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ABSTRACT 

A  general  method  is  presented  for  constructing  a  location  estimator  which 
is  asymptotically  efficient  at  any  two  different  location-scale  families  of 
symmetric  distributions  as  well  as  at  an  appropriately  defined  class  of 
distributions  lying  in  between.  The  method  works  by  embedding  the  two 
families  in  a  comprehensive  parametric  model  and  identifying  the  estimator 
with  the  MLE.  The  case  when  the  families  are  Normal  and  Double  exponential  is 
examined  in  detail. 
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SIGNIFICANCE  AND  EXPLANATION 


This  paper  considers  the  following  probl  in  the  estimation  of  the 
center  of  a  aymnetric  probability  distribution.  Suppose  the  statistician  has 
a  model  F  which  he  hopes  is  a  good  approximation  to  the  true  underlying 
distribution  H  generating  his  data.  Further  suppose  that  he  has  reason  to 
believe  that  any  deviation  of  H  from  F  will  probably  be  in  the  direction 
of  another  model  G.  A  general  procedure  is  presented  for  constructing  an 
estimator  which  is  asymptotically  efficient  at  both  F  and  G  as  well  as  at 
a  suitably  defined  family  of  distributions  lying  in  between.  The  case  when 
F  is  Normal  and  G  Double  exponential  is  studied  in  detail  via  both 
asymptotic  theory  and  Monte  Carlo  simulation  (for  finite  sample  sizes).  The 
estimator  is  shorn  to  compare  favorably  against  nine  other  well  known 
conpetitors .  Computer  programs  are  included. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


DIRECTIONALLY  EFFICIENT  ROBUST  ESTIMATORS 
OF  LOCATION  VIA  EXPONENTIAL  EMBEDDING 


Wei-Yin  Loh 

1 .  Introduction 

In  robust  estimation  of  the  center  of  a  symmetric  distribution,  we 
usually  assume  that  we  have  a  parametric  model  F  (e.g.  Normal  location-scale 
family)  which  we  hope  is  a  good  approximation  to  the  true  underlying 
distribution,  but  we  do  not  assume  it  to  be  exactly  right.  A  robust  estimator 
is  then  desired,  i.e.  one  which  is  efficient,  or  nearly  efficient,  at  F  and 
has  reasonably  good  efficiency  in  a  neighborhood  of  F.  In  the  case  that  F 
is  Normal,  the  neighborhood  is  guite  often  taken  to  consist  of  all  symmetric 
distributions  with  tails  ranging  in  thickness  from  the  Normal  to  the  Cauchy. 
Typically  no  particular  distribution  in  the  neighborhood  (other  than  F)  is 
preferred  over  the  others,  i.e.  we  do  not  require  the  estimator  to  be  more 
efficient  at  some  distributions  than  at  others. 

In  this  paper,  the  case  is  considered  where  the  statistician  has  reason 
to  believe  that,  if  the  true  distribution  were  to  deviate  from  F,  it  would 
probably  (but  not  definitely)  be  towards  a  heavier-tailed  model  G.  In  such 
circumstances  it  would  be  desirable  for  the  estimator  to  possess  high,  if  not 
optimal,  efficiency  at  G  as  well.  Gastwirth  (1966)  and  Crow  and  Siddiqui 
(1967)  have  studied  this  problem  when  G  consists  of  one  or  more  parametric 
families.  Both  these  papers  suppose  it  is  known  that  the  population  sampled 
belongs  to  a  set  F  of  parametric  families,  like  {Normal,  Double 
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exponential}  or  {Normal,  Double  exponential,  Cauchy).  They  then  search  for  an 
estimator,  within  various  classes  (like  R-  and  L-  estimators),  which  has 
maximin  efficiency  over  F.  Since  maximin  estimators  emphasize  safety  over 
efficiency  however,  they  may  be  efficient  (or  even  asymptotically  efficient) 
nowhere  in  F. 

If  asymptotic  efficiency  were  the  only  criterion,  of  course  solutions  to 
our  problem  are  available  from  the  class  of  fully  adaptive  estimators  (see 
e.g.  Andrews  et  al.  (1972),  Sacks  (1975),  Stone  (1975)  and  Beran  (1978)), 
these  being  constructed  to  be  asymptotically  efficient  at  all  sufficiently 
smooth  distributions.  But  since  these  estimators  make  no  use  of  our  knowledge 
of  F  and  G,  it  seems  plausible  that  a  semi-adaptive  procedure  which  uses 
this  information  explicitly  may  perform  better  (at  F  and  G)  for  small 
samples.  Examples  of  such  procedures  have  been  suggested  by  Hogg  and 
others.  Hogg  (1967)  proposes  an  estimator  which  chooses  among  the  sample 
mean,  median,  25%  trimmed  mean,  and  mid-range  according  to  the  sample 
kurtosls.  This  estimator  is  asymptotically  efficient  at  the  Normal,  Double 
exponential  and  uniform  distributions,  and  empirical  evidence  (see  e.g. 

Andrews  et  al.  (1972)  and  Wegman  and  Carroll  (1977))  suggests  that  its  small 
sample  performance  at  these  distributions  is  better  than  many  of  the  blatantly 
adaptive  procedures.  Further  it  is  robust  since  it  yields  the  sample  median 
or  trimmed  mean  whenever  the  sample  kurtosis  is  large.  Another  procedure  that 
applies  maximum  likelihood  or  Bayes  rule  to  first  select  a  family  of 
distributions  from  within  a  predetermined  set  and  then  uses  the  optimal 
estimator  for  that  family  is  considered  in  Hogg  et  al.  (1972). 

In  this  paper  we  introduce  a  general  method  for  constructing  a  semi- 
adaptive  estimator  which  is  asymptotically  efficient  at  any  two  parametric 
families  F  and  G  as  well  as  at  a  family  of  distributions  lying  between 


them.  The  definition  of  this  family  and  an  explanation  for  its  relevance  are 
given  in  Section  2.  Unlike  the  Hogg-type  procedures  which  are  based  on 
trimmed  means,  our  estimator  is  more  like  an  M-estimator  since  it  is  the 
maximum  likelihood  estimator  (MLE)  corresponding  to  a  genuine  likelihood 
function.  If  the  tails  of  F  and  G  span  a  wide  enough  range,  we  expect 
this  estimator  to  be  robust.  The  particular  example  when  F  is  Normal  and 
G  Double  exponential  is  studied  in  detail  in  Sections  3-5.  Section  3  derives 
the  likelihood  equations,  section  4  deals  with  the  asymptotic  properties  of 
the  estimator,  and  section  5  contains  empirical  evidence  on  its  small  sample 
performance  compared  to  some  well-known  adaptive  and  nonadaptive  competitors. 


2.  Some  embedding  methods 

Let  F  *  {0  ^(0  %x-0))j  -<•  <  0  <  •,  0  >  0}  and  G  =  {t  ^gft  '(x-©)); 

-«  <  0  <  •,  T  >  0}  be  two  location-scale  families  of  densities  on  the  real 
line.  We  assume  that  f(x)  and  g(x)  are  symmetric  about  0  and  want  an 

A 

estimator  0  of  0  that  is  asymptotically  efficient  at  F  and  G  as  well 
as  those  densities  in  between.  There  appears  to  be  no  universal  agreement  on 
what  is  meant  by  the  set  of  distributions  "between"  two  families.  We  will 
define  it  in  the  following  way.  Let  H  be  a  comprehensive  parametric  model 
parametrized  by  an  additional  parameter  X  e  [X^X^)  such  that  X  =  X^ 
corresponds  to  F  and  X  =  X^  corresponds  to  G.  Then  all  distributions 
in  H  that  are  not  in  F  or  G  will  be  considered  as  being  in-between 
and  G.  In  order  for  the  estimation  problem  to  be  well-defined,  the  densities 
in  H  will  have  to  be  symmetric.  Once  a  suitable  embedding  H  is  found,  we 
will  define  0  as  the  MLE  over  H  of  the  center  of  symmetry  0. 

There  are  at  least  three  approaches  to  constructing  such  an  embedding. 

The  linear  embedding  consists  of  densities  defined  by 


(2.1) 


h(x;0,O,T,X)  =  o  %1-X)f{o  *(x-0)}  +  T  *Xg{x  (x-0)}; 

-«  <  0  <  •»  0 , T  >  0;  0  <  X  <  1 


This  construction  has  the  advantage  of  being  simple  and  allowing  a  physical 
interpretation  for  X  as  a  prior  probability.  Unfortunately,  since  a  and 
x  are  unknown  parameters,  the  likelihood  function  is  unbounded  at  each  data 
point  Xj,«..,xn  for  n  >  2.  Hence  MLEs  do  not  exist.  This  difficulty  can 
be  avoided  by  including  the  restriction  that  o/t  =  constant.  But  besides 
drastically  diminishing  the  size  of  the  embedding,  this  seems  to  make  the 
choice  of  the  constant  artificially  important. 

Another  approach  is  to  have  the  members  of  H  be  "F  in  the  middle 


and  G  in  the  tails",  i.e.  for  k  >  0  define 


-4- 


a  ^{o  1  ( x— 0  )> 


I  x-0 I  <  k 


h(x;  0,0, T,k) 


T  ^{t  ^x-©)}  ,  otherwise 


The  M-estimator  of  Huber  (1964)  used  this  construction  with  F  Normal  and  G 
Double  exponential,  and  the  additional  requirement  that  h  be  continuously 
differentiable  in  x.  The  problem  of  estimating  all  the  parameters 
(0,0,T,k)  via  maximum  likelihood  does  not  appear  to  have  been  attempted, 
although  Bell  (1980)  has  investigated  the  question  of  adaptively  choosing  k 
from  the  data  using  the  criterion  of  minimum  estimated  asymptotic  variance. 

A  third  construction  is  the  exponential  embedding  where 

h(x;  0,0, x,X)  =  c(0,T,X)f1  ^{o  \x-0)}g*{T  ^x-O)}  ; 

(2.2) 


-•  <  0  <  •;  9,1  >  0;  0  <  )  <  1  ; 

and  c(0,T,X)  is  a  scaling  factor.  Cox  (1961),  Atkinson  (1970),  Brown  (1971) 
and  Weerahandi  and  Zidek  (1978)  have  used  this  in  various  contexts.  Recently, 
using  the  Kullback-Leibler  information  number  as  a  measure  of  statistical 
distance,  Loh  (1983)  showed  that  as  X  ranges  from  0  to  1,  the  distri¬ 
butions  represented  in  (2.2)  in  fact  constitute  the  shortest  path  between  f 
and  g  in  distribution-space.  This  result  offers  a  justification  for  claim¬ 
ing  that  (2.2)  yields  densities  in  between  F  and  G. 

We  adopt  the  method  of  exponential  embedding  in  this  paper  and  estimate 

A 

0  with  its  MLE  0,  regarding  o,  T  and  X  as  nuisance  parameters.  (The 
fact  that  the  nuisance  parameters  may  not  all  be  identifiable  is  not  worrisome 

A 

because  we  are  only  interested  in  estimating  0.)  It  is  clear  that  0  is 
location  and  scale  equivariant,  i.e.  if  we  transform  the  data  vector  )t  to 
aj<  +  b  for  some  constants  a  and  b,  then 

A  A 

0(ax+b)  =  a0(x)  +  b 

A  A 

If  0_  and  0„  are  the  MLFs  for  0  under  the  submodels  F  and  G, 

F  G 

A 

the  following  theorem  which  will  he  used  later  gives  conditions  for  8  to  lie 

A  A 

between  0_  and  0  . 

F  G 


-5- 


Theorem  2.1.  Suppose  that  for  each  (<j,t)  the  likelihood  functions  (of  0) 

A 

under  F  and  G  are  unimodal.  Then  the  MLB  0  for  (2.2),  when  0,  T  and 

A  A 

X  are  treated  as  nuisance  parameters,  always  lies  between  0  and  0  . 

F  G 

Proof.  This  follows  from  the  assumption  of  unimodality  of  the  likelihoods  and 

A  A 

the  fact  that  the  values  of  0_  and  0„  are  unchanged  whether  o  and  t 

F  G 

are  known  or  not. 

When  F  is  Normal  and  G  Double  exponential,  the  MLEs  are  the  sample 

mean  X  and  median  X  respectively.  These  two  estimators  are  at  the 

extremes  of  the  family  of  symmetrically  trimmed  means,  and  the  a-trimmed  mean 

X  is  often  thought  of  as  a  compromise  if  the  true  underlying  distribution  is 
a 

believed  to  have  tails  between  those  of  the  Normal  and  Double  exponential.  It 

A 

is  interesting  to  note  that  X^  does  not  share  the  property  of  0  in  Theorem 
2.1  in  this  case. 

A 

Since  0  is  an  MLB,  classical  theory  suggests  that  under  regularity 
conditions,  it  is  asymptotically  efficient  when  0  <  X  <  1.  When  X  =  0  or 
1,  proofs  of  asymptotic  normality  are  more  difficult  since  the  true  parameter 
vector  is  now  a  boundary  point.  There  appears  to  be  no  general  theorems  for 
such  situations.  In  specific  cases,  a  proof  will  probably  depend  on  a 
combination  of  the  results  of  Huber  (1967)  and  ad  hoc  arguments.  Robustness 

A  A 

of  0  is  likely  to  depend  on  the  robustness  of  0  if  6  is  heavier-tailed 
than  F.  Intuitively  this  is  because  the  estimated  density  (2.2)  will  tend  to 
be  close  to  some  member  in  G  if  the  true  distribution  is  heavier-tailed 
than  G.  we  illustrate  these  points  with  an  example  in  the  following 


sections. 


w 


>.  Normal-Double  exponential  example:  likelihood  equations 

We  study  here  the  case  when  F  is  Normal  and  6  Double  exponential. 
Equation  (2.2)  becomes 

12  2 

(3.1)  h(x»  0,s,t)  *  c(s,t)exp{-  —  s  (x-0)  -  t|x-0|) 


r 


where  c(s,t)  ■ 

j  s*(t/s)/*(-t/s) 

if 

s,t 

>  0 

< 

s/(2it)1/2 

if 

t  = 

0,  s  >  0 

t/2 

if 

s  = 

0,  t  >  0 

and  $(  ) ,  ♦(  )  are  the  standard  Normal  density  and  cumulative  distribution 
functions,  the  scale  factor  c(s,t)  is  defined  on  V  = 

[0,-)  x  [0,»)  \  {(0,0)}.  It  can  be  checked  (using  e.g.  (3.6)  below)  that 
c(s,t)  is  continuous  on  V.  Clearly  (3.1)  yields  F  when  t  =*  0  and  G 

when  s  •  0. 

AAA 

Let  0,  s,  t  be  the  MLEs  for  0,  s,  t.  the  following  theorem  whose 
proof  is  sketched  in  the  Appendix  shows  that  the  three-parameter  minimization 
problem  of  determining  the  MLEs  may  be  reduced  to  one  involving  only  one 
parameter  in  [0,«). 

theorem  3. 1 .  Let  (Xf,#..*xn)  be  an  ordered  sample  of  size  n  and  x, 

*  “  (*,,  +  xr  /i-m*  *>e  the  sample  mean  and  median  respectively  (here 

lln+i  )/zj  in/z+ij 

[  ]  is  the  greatest  integer  function).  Suppose  that  x  <  x.  If  n  is  even 

and  xn/2  *  x»  then  0  *  x.  Otherwise  the  minimization  of  the  likelihood 

corresponding  to  (3.1)  can  be  reduced  to  the  following  one-dimensional 

2 

problem:  Let  v  ■  t/s,  w  *  t/s  and  let  kQ  be  the  largest  integer  k  such 

that  x.  <  x.  Define  to  be  .e  set  of  integers  (k  ,...,[ (n-1 )/2] }  and 
k  0 

divide  the  interval  [0,*j  into  the  subintervals  Up1k#P2k)»  ^qlk'q2k^f 
k  e  N}  where  p1k  »  0,  <12,  [(n-1  )/2]  =  ®  and  for  k  "  k0  +  •  •  •  *  l (n-3)/2) , 


-7- 


p_,  «  n(n-2k)  ’(x  -x){(n-2k)  '  (x  -x)I|x  -x 

2k  k+1  k+1  i  k+1 


(3.2) 


-1_,  ,2. -1/2 
+  n  I  (x.  -x  )  } 
i  k+1 


q1k  *  p2k  * 

q2k  =  n(n-2k“2>  1 (xk+1~x){ (n-2k-2)  1  (xk+1_x)E !xi~xk+1 


,  2.  -1/2 

+  n  r<xi-xk+1)  } 

p1,k+1  =  q2k  * 

For  each  k  e  N,  define  for  v  e  [0,«)  the  function 

f  2  2*1^  # 

(3.3)  fk(v)  ~  j  log(w/v)  +  log  •(-v)  +  v  (2nw  )  {][  (x^O-w)* 


I  (x.-8-Kf)  } ,  if  ve  (p  ,q,  J 
k+1 


0  ,  otherwise 


where 


(i)  for  ve  tp1k#p2k^» 

-1  2  £  -  ,  ^  ,  -2  4.r  -  ..2  -1  2_ .  -.2.1/2 

w  =  n  v  )  (x-x.)  +  {n  v  {£  (x-x  )}  +  n  v  E(x  -x)  } 

1  i  1  l 


(3.4)  0  =  x  +  w(n-2k)/n  ; 

(ii)  for  v  e  t<J1k*q2kJ » 


w  «  (2n)  1v2E|x.-x  |  +  n  V(I|x  *x  J)2 

l  k+1  2  i  k+1 


x  2r/  .2,  V2 

+  4n  v  Hxi-xk+1)  } 


(3.5)  0  = 


Vi  • 


r>  w  _ 

Let  f(v)  =  l  f  (v)  and  v  minimize  f(v)  over  [0,«).  Then  the  0 
keN  . 

corresponding  to  v  (given  by  (3.4)  or  (  ’.5))  is  the  VLE  for  0  for  the 


density  (3.1) 


We  are  unable  to  prove  or  disprove  that  0  is  a.s.  unique,  although  our 


experience  with  numerical  examples  suggests  that  this  is  the  case.  Should 

multiple  roots  occur,  however,  we  can  choose  the  root  closest  to  the  sample 

median.  In  view  of  the  theorems  in  the  next  section,  this  will  guarantee  a 

consistent  sequence  of  roots.  To  implement  the  method  on  a  computer,  we  note 

that  if  v*  lies  in  [q,  ,,  its  value  need  not  be  computed 

i ,  l  ( n-  i )  / 1 J 

A 

exactly  since  0  is  independent  of  it.  So  we  need  only  determine  whether 

* 


v  lies  in  this  interval.  This  search  is  greatly  assisted  by  the  following 

approximation  which  effectively  reduces  the  interval  to  a  finite  one. 

2 

Theorem  3.2.  Let  S,  =  E|x.-x.  .1  and  S„  =  E(x,-x,  . „ )  where  k  = 
-  1  i  k+1  2  l  k+1 

[  ( n- 1 )  /2  ]  .  If  v  >  max{2/2,  2S^VnS^,  q^},  then  f(v)  -  log(  S  y/(  n/iiT)  )  -  1  - 

-2  2  -4  2  2 

v  {nS^/^S*)  -  1}  is  bounded  above  by  v  {5/2  +  (7/12 )  (nS2/S1  )  }  and  below 

by  -v“4{(3/2)(nS2/S^)2  -  (*/2)(1  -w/8)}. 

_2 

Proof .  Expand  f(v)  in  powers  of  v  using  Taylor's  series  and  the 
inequalities  (see  e.g.  Johnson  and  Kotz  (1970)  p.  279):  for  x  >  0, 

(3.6)  {(x2+8)V2  +  3x}/4  <  <|>(x)/*(-x)  <  {(x2+2*)1/2  +  (n-1)x}/ir  . 

A  Fortran  program  that  uses  these  two  theorems  is  given  in  the 
Appendix.  It  uses  a  modified  version  of  the  function  minimization  routine 
FMIN  in  Forsythe,  Malcolm  and  Moler  (1977)  and  calls  IMSL  subroutines  MDNORD 
and  MSMRAT  to  compute  $(x)  and  Mill's  ratio. 


It  is  shown  in  this  section  that  8  is  asymptotically  efficient  when  the 
model  (3.1)  is  correct.  If  (8,s,t)  is  an  interior  point  of  the  parameter 
space,  standard  methods  can  be  used  to  prove  consistency  and  asymptotic 

AAA 

efficiency  of  (0,s,t).  However  these  methods  are  inapplicable  when  (0,s,t) 
is  a  boundary  point,  as  when  the  true  underlying  distribution  is  Normal  or 
Double  exponential.  For  these  cases  we  use  a  thoerem  in  Huber  (1967)  to  prove 
consistency  and  then  resort  to  ad  hoc  methods  to  argue  asymptotic  efficiency. 
Incidentally,  consistency  of  0  alone  is  a  consequence  of  Theorem  2.1  since 
it  is  sandwiched  between  X  and  X  both  of  which  are  consistent  under 

A 

(3.1) .  Consistency  for  the  other  MLEs  (as  well  as  9)  is  shown  in  the 
following  theorem. 

A  A  A 

Theorem  4.1.  0 ,  s  and  t  are  consistent  estimates  of  0 ,  s  and  t  when 

(3.1)  is  correct. 

Proof .  The  proof  consists  of  checking  that  the  conditions  in  Theorem  1  in 
Huber  (1967)  holds.  These!  are  called  assumptions  (Al)  -  (A5)  in  the  paper  to 
which  we  refer  the  reader  for  a  precise  statement.  Let  ft  =  i—0,*0)  x  V, 
where  D  is  defined  in  Section  3.  Let  o  =  (8,s,t)  e  ft,  and  = 

(0o,Vto)  be  the  parameter  vector  corresponding  to  the  true  distribution. 
Following  the  suggestion  in  Huber  (1981,  p.  130),  we  take  pairs  yn  = 
(x2n_1,X2n>  of  the  original  data  (x^,x2,...)  as  our  new  observations  and 
define 

p(y,a)  =  -2  log  c(s,t)  +  ~  s2{(x1“0)2  +  (x2~0)2} 

+  t{ | x 1 -9  I  +  Ix2-0|}  . 

Assumption  (A-1)  is  immediate,  and  (A-2)  follows  from  the  continuity  of 

p(y,a)  as  a  function  of  a.  Let  a(y)  =  p(y,ag)  and  note  that 

E{p(y,a)  -  a(y)}  is  a  Kullback-Leibler  information  number;  hence  it  is  non- 


negative  and  well-defined  (possible  ■»•)  for  all  a  e  Jl,  and  vanishes  only 

when  a  *  a^.  (Here  all  expectations  are  taken  with  respect  to  a^.)  This 

implies  assumptions  (A-3)  and  (A-4). 

Finally  to  verify  (A-5),  let  «  be  the  point  at  infinity  in  the  one- 

point  compactification  of  fl.  In  our  context  (A-5)  may  be  stated  as  follows: 

There  is  a  continuous  function  b(a)  >  0  such  that 

(i)  inf{p(y,a)  -  a(y)}/b(a)  >  h(y)  for  some  integrable  h, 
a 

(ii)  lim  inf  b(a)  >  0,  and 

O-MJ 

(iii)  E{lim  inf  [p(y,a)  -  a(y)3/b(a)}  >  1. 

a-M» 

Take  b(a)  to  be  identically  1  for  all  a.  Then  (ii)  is  immediate  and  (i) 

will  follow  if  we  show  that  inf  p(y,a)  is  integrable.  For  each  s  and  t, 

a 

p(y,a)  is  minimized  when  8  =  /2.  Therefore  writing  z  =  Ix^-Xjl/^, 

we  see  that 

2  2 

(4.1)  inf  p(y,a)  *  inf  {-2  log  c(s,t)  +  s  z  +  2tz} 

a  8ft 

Now  suppose  (s,t)  is  an  interior  point,  and  make  the  change  of  variable  u  = 

t/s.  The  expression  in  parenthesis  on  the  RHS  of  (4.1)  can  be  rewritten  as 

2  2  2 

(4.2)  H(z,8,u)  ■  -2  log  s  +  u  +2  log  8(-u)  +  s  z  +  2suz 

2  1/2 

For  fixed  u  this  is  minimized  when  2sz  =  (u  +4)  -  u.  Substituting  this 

for  s  in  (4.2)  and  differentiating  with  respect  to  u  yields 

2  1/2 

dH/du  =  u  -  2$(u)/$(-u)  +  (u  +4) 

which  is  positive  for  all  u  (see  Birnbaum  (1942)).  We  therefore  conclude 

that  inf  H(z,s,u)  is  attained  at  u  =  0  and  s  =  z  \  Hence  inf  p(y,a)  = 

a 

2  log  z  +  constant,  which  is  clearly  integrable,  and  so  (i)  obtains.  To 

verify  (iii)  it  suffices  to  check  that  lim  inf  p(y,a)  =  ».  There  are  two 

a-H» 
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cases  to  consider,  namely  (a)  |0|,  s,  t  ♦  •»,  and  (b)  |0|  ♦  »,  s,t  ♦  0.  In 
both  cases  however,  we  have 

z'  2  2  -1 

p(y,a)  *  |  0(-  log  s  +  s  0  /2  +  t|0|)  if  ts  ♦  constant 

^  0{-  log  t  +  s202/2  +  t|0|  )  if  ts"1  -*■  •  . 

Obviously  lim  inf  p(y,a)  =  •  in  either  case.  Aus  (iii)  is  verified  and  the 

proof  is  ended. 

We  are  now  ready  to  deduce  asymptotic  normality  and  efficiency. 

A 

Theorem  4.2.  0  is  asymptotically  efficient  when  the  underlying 

distribution  F  has  a  density  given  by  (3.1). 

Proof .  When  the  true  parameter  vector  (0,s,t)  is  an  interior  point  of  the 
parameter  space,  the  asymptotic  efficiency  of  the  MLEs  follows  easily  from  the 
standard  theorems  (see  e.g.  Lehmann  (1983)).  We  therefore  only  prove  the 
result  when  either  s  or  t  is  zero. 

*  P  *  P 

(i)  F  is  Normal.  By  the  preceding  theorem,  t  ♦  0  and  s  ♦  s  for  some 
s  >  0.  It  is  clear  that  equations  (3.4)  and  (3.5)  together  define  a  monotone 
function  of  v  in  [0,«).  Therefore 

v  A  —  A  A  a  ^ 

/n  |0-x|  <  w|n-2k|//n  =  t|n-2k|/(s  /n) 
for  some  k  between  kQ  (defined  in  Theorem  3.1)  and  n/2.  Here  |n-2kg| 
is  the  difference  between  the  number  of  deviations  {x^-x}  with  positive 
signs  and  the  number  with  negative  sign.  Since  |n-2kQ|  =  Op(/n)  (David 
(1962)),  we  see  that  /n  |0-x|  <  t|n-2kg|/(s  /n)  +  0.  Hence  0  has  the  same 
asymptotic  variance  as  x  which  is  efficient. 

(ii)  F  is  Double  exponential.  We  assume  without  loss  of  generality  that 

_  AW 

x  <  x  and  again  use  the  notation  of  Theorem  3.1.  We  know  from  Theorem  4.1 

Ap  A  p  A  A  A  p 

that  s  ♦  0  and  t  ♦  t  >  0.  Therefore  v  =  t/s  ♦  •.  Let  kn(x)  he  the 

A 

largest  integer  k  such  that  q2)c  (defined  in  (3.2))  satisfies  q2)c  <  v. 
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Then  kQ  <  k„  <  n/2  and  n-2kQ  =  Op(/n)  (see  Brown  and  Kildea  (1979)).  So 

(n-2k  )//n  converges  in  probability  to  some  random  variable  Y.  We  will  show 
n 

that  Y  is  degenerate  at  0.  First  observe  that  /n|x^  “xl  *  /n|x-x|  = 
k 

1  ^  ^  “■  p 

Op ( 1 )  and  n  £  (x-x.)  =  (2n)  E|x.-x|  +  constant.  Next  note  that  qj^ 

1  1  1 

can  be  written  as 

1  -  i  _  k+1  _ 

q,.  =  n(n-2k-2)~  (x.  .  -x) /{ 2( n-2k-2 )"  (x.. ,-x)  £  (x-x.) 

^2k  k+1  k+1  ,  i 


-1 

+  n 


—2, 1/2 

Kx^-x)  } 


*  P  P 

Since  v  +  •,  this  yields  q^  +  •».  This  implies  that  Y  =  0.  Therefore 

.  —  p  n  *  *.  * 

k  /n  =  —  +  r  where  /n  R  +0.  Since  x.  <  0  <  x,  it  follows  that  9  has 

n  2  n  n  k 

__  n 

the  same  asymptotic  variance  as  x  (Lehmann  (1983),  Chapter  5,  Problem  3.5). 
The  next  theorem  shows  that  for  heavy-tailed  distributions  like  the 

A 

Cauchy  (or  Tukey's  "slash")  0  is  asymptotically  equivalent  to  the  median. 

It  is  therefore  robust  against  these  distributions. 

Theorem  4.3.  Assume  that  the  true  underlying  density  has  tails  of  order 

—2  *  ^ 

|x|  .  Then  0  has  the  same  asymptotic  variance  as  X. 

Proof .  Assume  as  before  that  the  x's  are  ordered  and  use  the  notations  of 
Theorem  3.1.  Brown  and  Tukey  (1946)  showed  that  for  such  distributions 
(4.3)  X  =  Op(1),  n"1  E|X.|  =  0p(1),  n"1  EX2  *  Op(n)  . 


It  is  easy  to  check  that 

s"2  =  Op(n_1  E(x.-0)2),  v/s  *  Op(n_1  E|xi~0|)  . 

A  A 

A* 

Since  0  lies  between  X  and  X,  (4.3)  implies  that  v  *  Op(/n).  Now 

choose  {k  }  such  that  k  ■  ~  n  +  nR  where  /n  R  +0.  This  ensures  that 
n  n  2  n  n 

the  kfc>1  order  statistic  X.  has  the  same  asymptotic  variance  as  X. 
n  Kn 

Putting  V  *  in  (3.2)  we  see  also  that 
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q„  .  -  (/n  R  )  (x.  -x)/{(/n  R  )  (x  -x)(n  E|x  -x  |) 

2,k  n  K  n  x  x  k 

n  n  n  n 


-2  .2.1/2 

+  n  E(xi-xfc  )  } 

n 


which  in  view  of  (4.3)  is  at  most  0^((/n  R  )”1).  Clearly,  we  can  choose 

P  n 

* 

k  so  that  /n  R  ♦  0  and  nR  ♦  •.  Then  for  large  n,  q,  v  «  v  with 
n  n  n  ,Kn 

a  ^ 

high  probability.  This  yields  X^  <  8  <  X  and  hence  the  result. 


Small-sample  behavior:  sensitivity  curve,  breakdown  bounds  and 


Figure  5.1.  Sensitivity  curve 


X 


3  sample  median 


Adaptive; 


JOH 


John's  adaptive  estimator 


TAK  “  Takeuchl's  adaptive  estimator 
JAE  *  Jaeckel's  adaptive  trimmed  mean 
HG1  ^ 

>  Hogg-type  adaptive  estimators. 

HG2  J 

The  definitions  of  JOH,  TAK  and  JAE  are  given  in  Andrews  et  al.  (1972)  under 

the  same  names.  HG1  was  first  suggested  by  Hogg  (1974)  and  is  defined  to  be 

X  if  Q  <  1.81,  X  __  if  1.81  <  Q  <  1.87,  and  X  if  Q  >  1.87. 

Here  Q  is  the  ratio  {U(.2)  -  L( .2)}/{U( .5)  -  L(.5)},  where  U(B)  and 

L(0)  are  the  averages  of  the  largest  and  smallest  t(n+1)0]  order  statistics 

respectively.  The  estimator  HG2  is  a  modification  of  HG1  to  make  it 

asymptotically  efficient  at  the  Normal  and  Double  exponential  distributions. 

It  is  defined  to  be  X  if  0  <  1.81,  X  if  1.81  <  Q  <  1.87,  and  X  if 

0  >  1.87.  From  Table  5.1  it  is  clear  that  for  n  *  10,  20  or  40,  the 

* 

breakdown  bounds  for  8  are  superior  to  all  the  others  except  those  for  the 
sample  median. 

Monte  Carlo  estimates  of  the  variances  (multiplied  by  n)  of  each  of 
these  estimators  are  given  in  Tables  5.2  -  5.4  for  n  =  10,  20  and  40  and 
the  following  distributions:  (i)  N(0,1)  (Normal  with  mean  0  and  variance 
1),  (ii)  the  density  (2.2)  with  f  *  N(0,1),  g  =  Double  exponential  with 
density  j  e  ^ ,  and  A  =  ~,  (iii)  Double  exponential,  (iv)  contaminated 
Normal:  90%  N(0,1)  +  10%  N(0,100),  and  (v)  Cauchy.  The  first  three 
distributions  are  picked  for  the  study  because  they  span  the  range  in  which 

A 

8  is  efficient.  Distributions  (iv)  and  (v)  are  included  to  test  its 
robustness  properties.  Estimates  of  the  standard  errors  are  given  in 
parentheses,  and  the  minimum  estimated  variance  for  each  distribution  is 
underlined.  The  simulations  were  done  on  a  VAX/11/750  computer.  The  IMSL 
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generator  GGUtf  was  used  to  generate  uniform  random  numbers  and  the  Box-Muller 
transformation  applied  to  produce  normal  deviates.  The  Princeton  swindle  was 
used  whenever  possible.  The  number  of  replications  ranged  from  1000  -  5000. 

It  is  immediately  clear  from  these  tables  that  in  none  of  the  sample 

* 

size-distribution  combinations  considered  does  6  beat  all  of  the  other 
estimators.  To  analyse  them  further,  we  can  look  at  deficiencies.  These  are 
defined  to  be  the  ratios  (estimated  variance) /(minimum  estimated  variance), 
where  the  minimum  is  taken  over  the  ten  estimators  compared.  The  coded 
deficiencies  are  shown  in  Tables  5.5  -  5.7.  Only  deficiencies  greater  than 
1.5  appear  as  digits  or  x's.  All  the  estimators  (with  the  possible  exception 
of  the  mean  and  median)  seem  to  be  equally  good  at  the  Normal  and  Double 
exponential  distributions.  The  excellent  performance  of  TAK  at  distribution 
(ii)  also  stands  out. 

Finally  to  conqpare  the  relative  performance  of  the  estimators  over  all 
situations  combined  we  follow  Tukey's  (1979)  suggestion  to  look  at  maxima  and 
sums  of  deficiencies.  For  each  sample  size,  let  A(i)  and  B(i)  be  the  maximum 
and  total  deficiency  of  the  ith  estimator  over  a  set  of  distributions.  The 
estimators  are  then  ranked  according  to  the  values  of  {A(i)}  and  (B(i)}. 
These  two  criteria  are  denoted  by  "minimax''  and  "total"  in  Tables  5.8  -  5.9 
where  the  estimators  are  ranked  first  for  distributions  (i)  -  (iii),  and  then 
again  for  all  five  distributions. 

The  following  points  may  be  made  from  these  two  tables: 

(a)  For  n  ■  20  or  40,  TAK  and  JOH  appear  hard  to  beat.  For  n  =  10 
however,  the  picture  is  quite  different.  Here  JOH  is  somewhat  below  average 
when  only  distributions  (i)  -  (iii)  are  considered,  and  TAK  has  a  poor  showing 
for  all  distributions  combined.  The  reason  may  be  that  these  two  estimators 
are  over-adapting  at  this  sample  size. 


(b)  There  seems  to  be  little  to  choose  between  6  and  HG2.  Both  are 
consistently  good  over  all  three  sample  sizes . 

(c)  HG2  is  superior  to  HG1  for  the  distributions  considered. 

A 

(d)  The  adaptive  trimmed  mean  JAE  trails  6  and  HG2  almost  every  time. 

(e)  None  of  the  nonadaptive  trimmed  means  (including  the  mean  and  median)  are 
competitive. 

The  above  results  encourage  us  to  feel  that  for  the  kind  of  situation 
described  in  the  introduction,  our  proposed  procedure  will  produce  viable 
estimators. 
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6 

1.069 

0.416 

1.59 

1.98 

5.09 

(.004) 

(.009) 

(.02) 

(.03) 

(.16) 

X 

1.000 

0.347 

1.99 

10.17 

768180.0 

(.007) 

(.04) 

(.40) 

(73179.0) 

\10 

1.051 

0.460 

1.61 

2.78 

14.63 

(.002) 

(.009) 

(.03) 

(.13) 

(3.51) 

*.25 

1.157 

0.630 

1.41 

1.57 

4.15 

(.006) 

(.011) 

(.03) 

(.02) 

(.15) 

A* 

X 

1.382 

0.874 

1.46 

1.74 

3.36 

(.015) 

(.015) 

(.02) 

(.20) 

(.08) 

JOH 

1.194 

0.445 

1.62 

1.82 

4.38 

(.007) 

(.010) 

(.03) 

(.19) 

(.23) 

TAX 

1.048 

0.307 

1.71 

2.62 

18.55 

(.002) 

(.007) 

(.03) 

(.19) 

(5.42) 

JAE 

1.081 

0.428 

1.56 

1.87 

6.53 

( .004) 

(.009) 

(.02) 

(.13) 

(.53) 

HG1 

1.119 

0.513 

1.48 

1.67 

4.18 

( .005) 


(.010) 


(.02) 


(.06) 


(.18) 


Table  Variance  x(n=^oj 


N(0,1) 

X"  2 

DEXP 

1 0  %  1  ON 

Cauchy 

* 

6 

1.074 

0.373 

1.46 

1.97 

3.41 

(.005) 

(.008) 

(.02) 

(.02) 

(.09) 

X 

1.000 

0.349 

1.92 

11.40 

19921.0 

(.007) 

(.04) 

(.41) 

(14362.0) 

*,10 

1.056 

0.476 

1.52 

2.08 

7.81 

(.002) 

(.009) 

(.03) 

(.08) 

(.40) 

*.25 

1.190 

0.689 

1.30 

1.52 

3.27 

(.008) 

(.013) 

(.03) 

(.01) 

(-14) 

A* 

X 

1.494 

1.012 

1.31 

1.83 

2.79 

(.028) 

(.018) 

(.02) 

(.02) 

(.06) 

JOH 

1.127 

0.300 

1.43 

1.44 

2.88 

(.005) 

(.008) 

(.02) 

(.01) 

(.09) 

TAX 

1.048 

0.218 

1.55 

1.42 

3.70 

(.002) 

(.006) 

(.02) 

(.02) 

(.16) 

JAE 

1.102 

0.404 

1.41 

1.48 

3.58 

( .005) 


(.009) 


(.02) 


(.01) 


(.12) 


1.000 


0.350 

(.007) 


31 

02 


1.99 

(.05) 


.5 


,95 

2 

,02) 

( 

49 

235438 

,33) 

(118967, 

,59 

6 

,03) 

( 

1.199 

(.009) 


1.513 

(.023) 


1.097 

(.004) 


1.035 

(.002) 


0.727 

(.014) 


1.166 

(.021) 


0.217 

(.005) 


0.139 

(.004) 


1.25 

(.03) 


1.23 

(.02) 


1.27 

(.02) 


1.42 

(.02) 


1.50 

(.01) 


1.87 

(.02) 


1.43 

(.01) 


1.35 

(.01) 


2.89 

(.06) 


2.62 

(.04) 


2.50 

(.06) 


2.73 

(.07) 


Table  S.S.  Deficiencies  of  estimators 


( n- 1 0 ) 


Variances  divided  by  minimum  variance  among  10  estimators  rounded 
to  nearest  integer.  One's  are  suppressed,  numbers  >  9  are  coded  x. 


Normal 


DEXP 


10% ION  Cauchy 


0 


2 


X 


X 


JOH 

TAX 

JAE 

HG1 

HG2 


2 

3 


2 


6 

2 


2 


x 

4 

2 

3 

e 

6 

2 

2 
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Table  5.7.  Deficiencies  of  estimators  ( 


n-40) 


Normal 


DEXP 


10% ION 


Cauchy 


Table  5.8.  Rank a  of  estimators  for  diets.  (1)  -  (iii ) 


n  Criterion  1 


3  4 


10  Minimax  TAX  HG2  6  JAE  X  JOH  X 


HG1  X 


Total  TAX  0  X  HG2  JAE  X  JOH  HG1  X  „ 

.10  .25 


Appendix 

Proof  of  Theorem  3.1.  Let  (xj,...,xn)  be  an  ordered  sample  of  size  n  such 
—  ~  2 

that  x  <  x.  Let  v  ”  t/s  and  w  »  t/s  .  Then  (3.1)  can  be  rewritten  as 

h(x;  0,v,w)  -  {2/2*  wv  *0(-v))  ^expt-  ~  v2w  2(|x-0|  +  w)2}  , 

-  •  <  0  <  •  ,  v,w  >  0  . 

First  observe  that  if  n  is  even  and  xn/2  <  x*  then  0  =  x.  So  for  n 

even,  we  need  only  consider  the  case  when  x  <  x  .  For  fixed  v  and  w, 

n/z 

the  likelihood  is  maximised  by  that  integer  k  and  0  satisfying 


(A.1) 


<  0  <  x. 


which  minimizes  J  (x. -0)  +  2w  J  fx  -0|.  It  is  clear  from  the  assumptions 

1  1  1  1 

have  made  and  Theorem  2.1  that  2k  <  n.  This  gives 

0  •  (  x  +  w(n-2k)/n  if  w  e  I. 


(A. 2) 


C  x  +  w( i 

V  XV*1 


if  w  e  J, 


where  1^  is  the  interval  ((n(xfc“x)/(n-2k)}  ,  n(xk+1-x)/(n-2k)]  and 


n( x.  .  -x) /<n-2k-2 ) ]  if  n  -  2k  -  2  >  0 
k+1  k+1 

tn(x.  -x)/(n-2k)  ,•)  if  n  -  2k  -  2  <  0 
k+1 


Now  for  fixed  v,  and  (0,k)  given  by  (A. 2),  it  can  be  verified  that  there 
is  a  unique  w  that  minimizes  the  likelihood.  This  w  is  given  by 


n  V  £(x-x  )  +  {n”Zv*(£(x-x  ))*  +  n  ’v2  £(x  -x)Z}V*  if  w  e  I 
1  1  1  1 


-  _  xx2  ^  -12  r,_  “»2,  1/2 


-12  ?, _  .  .  1  ,-24,?. _  i  »2  ■  ,-12  ? _  .2.  1/2 


(2n)"  v  5!  ^xi-x]c-+1  *  +  2  *n”  V  <^xi“xk+1*)  +  4n  v  ^(xi"xk+1)  * 


if  w  e  Jk  . 
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Inverting  these  equations  yields  the  intervals  ^  '  [q. *q.>:  ;  ; 

k  -  k0,...,[(n-1)/2]>  given  in  (3.2). 
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Fortran  programs 

SUBROUTINE  TSUB < XDATA , X , N # TOL , THETA ) 

C  XDATA  IS  AN  N-D IMENS IONAL  VECTOR  CONTAINING  THE  ORDERED  X-VALUES 
C  X  IS  AN  N-DIMENS IONAL  WORK  VECTOR 
C  TOL  IS  THE  PRECISION  REQUIRED  IN  FM IN 
C  THETA  CONTAINS  THE  COMPUIED  MLE 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 

DOUBLE  PRECISION  XDATA < N > . X ( N ) , VL < 2 f 50 ) , VR ( 2 , 50 ) 

DOUBLE  PRECISION  VO < 2 , 50 ) , SS ( 50 ) , SA ( 50 ) , SK < 50 ) , FV < 2 , 50 ) 

COMPLEX  ZSM , ZLG 

LOGICAL  GCASE 

COMMON  DTI ,RNf DRSyDRK 

EXTERNAL  F1,F2 

DATA  SQRT8/2.828427125D0/ 

RN~DBLE<  N ) 

NHALF=N/2 
NHP1-NHALF+1 
M*  <  N-l ) /2 
IFLAG*N-2ANHALF 
DO  5  1*1. N 
5  X  < I ) =XDATA  < I ) 

S=0.0 
T1 *0 . 0 
T2=0.0 

DO  10  1=1, NHALF 
T1=T1+X( NHP1 -  I ) 

T2=T2+X  <  NHALF + I ) 

.1.0  S=3+X(NHP1-  I)  AA2+X (NHALF*  DAA2 

IF < IFLAG  . EQ .  0)  THEN 
XMED=0.5A<X( NHALF) +X(NHP1 )  ) 

ELSE 

XMED=X ( NHP1  ) 

T2---T2+X ( N ) 

S~S+X  <  N ) AA2 
END  IF 

XBAR=  <  TI +  T2  > /RN 
AX2=S/RN 
TEMP=XBARAA2 
XVAR=AX2-TEMP 

SCEN=S-RNATEMP 
I F  <  X  B  A  R  .EQ.  XMED)  THEN 
THETA* XBAR 
GCASE*. FALSE. 

GO  TO  700 
END  IF 

IF ( XBAR  .LT.  XMED)  THEN 

GCASE*. FALSE. 

ELSE 

DO  20  1= J, NHALF 
TEMP*-*  <  I  ’ 
x  ( I  >  =-x  •'  n-  i>  i  > 

20  X(N-I+1  )*TEMP 

IF  ( IFLAG  .EQ.  M  X  <  NHP1  '  =  X  ( NHP1  > 

GCASE*. TRUE. 

XDAR=-XPAR 


-29- 


Apr  20  10:17  1983  t»xt  Pag»  2 


END  IF 
K0*NHALF 

30  IF(X(K0)  .LE.  XBAR)  GO  TO  35 

KO-KO-l 
GO  TO  30 

35  IF< ( IFLAG  .EQ.  0)  .AND.  <K0  .EQ.  NHALF ) >  THEN 

THETA-XBAR 
GO  TO  700 
END  IF 
81*0.0 

DO  50  I*KOf 1,-1 
50  S1*XBAR-X< I)+S1 

K*KO 

TEMP=0.0 
60  KP1-K+1 

0L<  1 ,K)*TEMP 
RM2K=DBLE(N-2AK> 

AS1 -SI /RM2K 

DIFF*X<KP1)-XBAR 

0NUN*DIFFARN/RM2K 

gK  sg  j 

0R<  1 , K)=0NUM/SQRT<2.0ADIFFAAS1+X0AR> 
UL(2,K)=0R<1,K) 

S-0.0 

DO  70  1=1 , K 

S*S+ABS< X <K+I>-X<KP1> ) +ABS<  X  <KP1- I>-X(KP1>  ) 

70  CONTINUE 

DO  80  I*2AK+1,N 

80  3=S+ABS<X( I)-X(KP1  ) ) 

SA  <  K  )  =S 

ASQ=AX2*X  <KP1  )A(X<KPl)-2.  OAXBAR  ) 

SS  <  K ) =ASQARN 

IF(K  .EQ.  M)  GO  TO  100 

RM2KP*DBLE<N-2AK-2> 

0NUM2*D IFFARN/RM2KP 

OR ( 2  f  K ) =0NUM2/SQRT  <  D  IFFAS/RM2KP+ASQ  > 
TEMP=0R<2,K) 

K*K+1 

S1=S1+XBAR-X  <K ) 

GO  TO  60 

100  DRS=RNASCEN 

DTEMP  =  0 . 5D0AL0G <  DBLE <  XOAR )  ) 

DO  210  K=KO,M 
DTI  =  SK  <  K ) 

DRK  =  DBLE<  K ) 

AA-OL  <  1 ,  K  ) 

BB*OR< 1 ,  K  > 

0X=FMIN<AA,BB,F1  .TOD 
00  < 1 ,  K )  =0X 

FV( 1 ,K )=F1 <  OX  >  «  DTEMP 
210  CONTINUE 

DTEMP2=L0G<5.0132565'18ARN) 

DO  220  K=KO,M-I 
DTI =SA ( K  > 

DRS=RNASS  <  K ) 

AA=0L(2.K> 
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BB»VR<2,K) 

VX*FMIN<AA.BBfF2.T0L> 

00<2fK)=UX 

FV(2fK)=F2<VX) -D TEMP 2 
220  CONTINUE 

FVM=FV  < 1 , KO ) 

Kl-KO 
11  =  1 

DO  230  K=K0, M 
DO  230  1  =  1,2 

IF<  FVH  .GE.  FV (  I , K ) )  THEN 
FVM«FV< I, K) 

K1*K 
11  =  1 
END  IF 

IF ( K  .EQ.  H)  GO  TO  235 
230  CONTINUE 

235  IF(K1  .EQ.  M  .AND.  V0<1, M>  .EQ.  VR(l.M))  GO 
DT1=SA<M> 

DRS=RNASS(M) 

AA=VR  <  1 ,  M ) 

BB=MAX  < SQRT8, 2 . OARNASQRT ( DBLE ( ASQ ) ) /DTI > 
00<2,M)=FMIN(AA,BB,F2,T0L) 

FV(2,M)=F2(V0<2,M) > -DTEMP2 
1=2 
K=M 

IF < FVM  .GE.  FV<2,M) )  GO  TO  300 
SRATI0=SSfM>ARN/SA(M)AA2 

C0N=L0G(SA<M)/2. 506628275 ARN ) +1 . 0 

A=CON-FVM 

B=SRATI0/2. 0-1.0 

C=0. 9539460517-1 .5A3RATIQAA2 

AA=BB 

IF(B >240,250 ,260 

!40  TEMP=5.0+7.0ASRATI0AA2/6.0 

VINF=SQRT(AHAX1 < -IEMP/B , 8 . 0 ) ) 
YU=CQN+R/VINFAA2+0.5ATEMP/VINFAA4 
IF ( FVM  .LE.  (CON+B/8.0+C/64.0) )  THEN 
GO  TO  600 

ELSE  IF  <  FUM  .GE.  YU)  THEN 

GO  TO  300 

ELSE 

CALL.  ZQADR(A,B,Cf Z8Mf ZLG, IER) 

IF (  IER  . NE  .  0)  NR  HE <  A,  A >  •  IER  FROM  ZUAIiR=’ 
BB=SQRT  <  MAX  <  DBLE ( ZSM ) ,  DBLE <  ZLG ) ) ) 

GO  TO  230 

END  IF 

150  IF<  A )  300  , 300 . 270 

160  IF(A)  300 , 265 ,270 

65  BB=-C/B 

GO  T'.  275 

!70  CALL  ZQADR <  A ,  B  , C .ZSM . ZLG . IER  > 

IF<  IER  .NE.  0)  WRITE(A,A)  •  IER  FROM  ZQAHR=* 
BB=SQRT<MAX(DBLF.(ZSM)  f DBLF. (ZLG)  )  ) 

IF ( BB  .LE.  AA )  GO  TO  600 
VX=FMIN(AA,BB.F2,T0L ) 


TO  300 


,  IER 


,  IER 

1 


75 

80 
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TEHP»F2 ( VX ) -DTENP2 

IF  <  F2  <  VX ) -DTEMP2  .GT.  FVM )  GO  TO  600 
300  THETA*X (H  +  l ) 

GO  TO  700 

600  IF( II  .EQ.  1)  THEN 

K=K  1 

g*vo<i,K> 

TEMP=VAA2/RN 

TEMP2*TEMPASK<K) 

U0*TEMP2+SQRT<TEMP2AA2+TEMPASCEN> 
TH£TA=XBAR+W0A<RN-DBLE<2AK) >/RN 
GO  TO  700 
ELSE 

THETA=X  <  K1  +  1 ) 

END  IF 

700  IF(GCASE)  THEN 
THETA=-THETA 
XBAR=-XBAR 
END  IF 
RETURN 
END 

DOUBLE  PRECISION  FUNCTION  FI ( V > 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 

COMMON  S1.RN.RS,RK 
U=VAS1 

D=U/SGRT  ( RS ) 

D2-DAA2 

ASSORT  <  1  . 0+D2 )  +I» 

CALL  MDNORD ( -V, P ) 

F1=LOG(A*P}+0.5/AAA2+2.0A<D/A+RKA<RN-RK)A<V/RN)AA2> 

RETURN 

END 

DOUBLE  PRECISION  FUNCTION  F2<V> 

IMPLICIT  DOUBLE  PRECISION  <A-H,0-Z> 

REAL  RM 

COMMON  SI ,RN ,RS 
A=4 . OARS 
U-VAS1 

Y- U+SQRT  <  UAA2+A ) 

CALL  M3MRAT ( REAL ( V  > , RM , IER ) 

IF <  IER  .NE.  0)  RM~  <  SORT <  VAA2+8 . DO  >  +3 . OAV ) /4 . DO 

F2-:L0G<  Y/RM  >  +2 .  OA  ( U+RS/Y )  /  Y 

RETURN 

PND 

DOUBLE  PRECISION  FUNCTION  FM  IN  <  AX  .  BX  .  F  .  TOI. ) 

C  THIS  IS  A  SLIGHTLY  MODIFIED  VERSION  OF  A  PROGRAM  BY  THE 
C  SAME  NAME  IN  FORSYTHE,  MALCOLM  AND  MOLER <3.3??  > 

DOUBLE  PRECISION  AX,BX,F.TOL 

DOUBLE  PRECISION  A ,  B  .  C  ..  D  .  E  ,  EPS  ,  XM  .  ?  ,  Q  ,  R  .  T0L1 ,  TCL2.U .  V  ,« 
DOUBLE  PRECISION  FU , FV , FU , FX , X 
SAVE  EPS 

D  A  T  A  C  /  0 . 3  8 1  9  6  6  0 1 1 3 D  0  / 

IF (EPS  .AT.  O.ODO)  GO  TO  15 
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EPS*1 .00 

10  EPS=EPS/2.00 

T0L1=1 .0+EPS 

IF ( T0L1  .GT.  1.00)  GO  TO  10 
EPS=SQRT  <  EPS ) 

15  A=AX 

B=BX 

V=A+CA ( B-A ) 

W=U 

x*v 

E=0.0 
FX~F<  X ) 

FV=FX 

FW=FX 

20  XM=0 . 5A ( A  +  B  > 

TOLl = EPS A ABS  <X)+T0L/3.0 
T0L2=2 . 0AT0L1 

IF ( ABO ( X-XM )  .LE.  <T0L2  -  0.5A<B-A>>>  GO  TO  90 
IF ( ABS ( E )  .LE.  TOLl)  GO  TO  40 
R*(X-W)A<FX-FV) 

Q=(X-V)A(FX-FW) 

P=(X-V)AQ-<X-W)AR 
0=2 . OOA  <  G-R ) 

IF  <  Q  .GT.  0.0)  P=-P 
Q=ABS ( Q ) 

R-E 

E=D 

30  IF < ABS < P )  .GE.  ABS ( 0 . 5AQAR ) )  GO  TO  40 

IF  <  P  .LE.  Q  A  <  A  -  X  ) )  GO  TO  40 
IF  <  P  .GE.  QA  <  B-X )  )  GO  TO  40 
D  =  P/0 
U=X+D 

IF ( <  U-A )  .LT.  T0L2)  D=S IGN( TOLl , XM-X ) 

I F  <  <  8 - U )  .LT.  T0L2 )  D=3  IGN ( TOLl , XM-X ) 

GO  TO  50 

40  IF ( X  .GE.  XM)  E=A-X 

IF ( X  .LT.  XM)  E~B-X 
D  =  CAE 

50  IF ( ABS ( D )  .GE.  TOLl)  U=X  +  D 

IF ( ABS  <  D )  .LT.  TOLl)  U  =  X  +  S IGN< TOLl ,  EO 

FU=F<U> 

I F  <  F  U  .GT.  FX  >  GG  TO  GO 

IF  <  U  .GE.  X)  A  -  X 

IF  ( U  .LT.  X)  B-X 

V  =  W 

FO-FW 

W  =  X 

FU-FX 

y=u 

FX~FU 
GO  TO  20 

d 0  IF  <  U  .LT.  A  ■  U 

IF ( U  .GE.  X)  B=U 
IF <  FU  .LE.  FW)  GO  TO  70 
1F(U  EQ .  X)  GO  TO  70 
IF '  FU  .  LE.  F'm  GO  TO  FO 
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IP(V  .E0.  X)  GO  TO  80 
IF<U  .EQ.  W>  GO  TO  80 
GO  TO  20 

70  U*W 

FV-FW 
U«U 
FW*FU 
GO  TO  20 

80  V=U 

FV*FU 
GO  TO  20 

90  T0L2*2.0AT0L 

IF< X-AX  .LT.  T0L2)  THEN 

IF<F(X>  .GE.  F <  AX  )  )  X=AX 

ELSE  IF  <BX-X  .LT.  T0L2)  THEN 

IF(F(X)  .GE.  F ( BX ) >  X»BX 

END  IF 

FMIN-X 

RETURN 

END 
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